clear all;
% clc;
warning off;

%% processing parameters
patval=1:5;
diamval=[3,6,9,12,15];
nopat=length(patval);
nodiam=length(diamval);
noaxons = 200;

isens_expsub = [0.20,0.35,0.70,0.55,0.80];
idisc_expsub = [0.50,0.50,1.50,1.20,3.10];
isens_expepi = [1.30,19.0,7.00,3.00,9.00];
idisc_expepi = [2.20,24.5,19.0,5.00,19.0]; 

%% Import Data 

mrg_epi_dcdata=cell(nopat,1);
mrg_sub_dcdata=cell(nopat,1);
sw_epi_dcdata=cell(nopat,1);
sw_sub_dcdata=cell(nopat,1);

iunit_epi_data=zeros(nopat,1);
iunit_sub_data=zeros(nopat,1);

% unit currents
for ii=1:nopat
    ptag=['p',num2str(patval(ii)),'_'];
    cd([ptag,'mrgpotentials']);
    iunit_epi_data(ii)=load([ptag,'epi_bp_current.txt']);
    iunit_sub_data(ii)=load([ptag,'sub_bp_current.txt']);
    cd ..;
end

% dc data
for ii=1:nopat
    
    % declare matrices to store incoming data
    mrg_epi_dcdata{ii}=zeros(noaxons,nodiam);
    mrg_sub_dcdata{ii}=zeros(noaxons,nodiam);
    sw_epi_dcdata{ii}=zeros(noaxons,nodiam);
    sw_sub_dcdata{ii}=zeros(noaxons,nodiam);
    
    ptag=['p',num2str(patval(ii)),'_'];
    
    % MRG data
    cd([ptag,'mrgthresholds\']);
    for jj=1:nodiam

        dtag = [num2str(diamval(jj)),'um'];
        tmp_a=load([ptag,'epi_bp_dc_act_',dtag,'.txt']);
        tmp_b=load([ptag,'sub_bp_dc_act_',dtag,'.txt']);
        mrg_epi_dcdata{ii}(:,jj)=abs(tmp_a(:,1))/(2/iunit_epi_data(ii)/1e3);
        mrg_sub_dcdata{ii}(:,jj)=abs(tmp_b(:,1))/(2/iunit_sub_data(ii)/1e3);        

    end
    cd ..;

    % SW data
    cd([ptag,'swthresholds\']);
    for jj=1:nodiam

        dtag = [num2str(diamval(jj)),'um'];
        tmp_c=load([ptag,'epi_bp_dc_act_',dtag,'.txt']);
        tmp_d=load([ptag,'sub_bp_dc_act_',dtag,'.txt']);
        sw_epi_dcdata{ii}(:,jj)=abs(tmp_c(:,1))/(2/iunit_epi_data(ii)/1e3);
        sw_sub_dcdata{ii}(:,jj)=abs(tmp_d(:,1))/(2/iunit_sub_data(ii)/1e3);        

    end
    cd ..;
    
end

clear tmp_a;
clear tmp_b;
clear tmp_c;
clear tmp_d;

%% plot ranges

yepi_min=1e6;
yepi_max=0;
ysub_min=1e6;
ysub_max=0;
for ii=1:nopat
    tmp_epi_min=min(min( [mrg_epi_dcdata{ii},sw_epi_dcdata{ii}] ));
    if(tmp_epi_min<yepi_min)
        yepi_min=tmp_epi_min;
    end
    tmp_epi_max=max(max( [mrg_epi_dcdata{ii},sw_epi_dcdata{ii}] ));
    if(tmp_epi_max>yepi_max)
        yepi_max=tmp_epi_max;
    end    
    tmp_sub_min=min(min( [mrg_sub_dcdata{ii},sw_sub_dcdata{ii}] ));
    if(tmp_sub_min<ysub_min)
        ysub_min=tmp_sub_min;
    end    
    tmp_sub_max=max(max( [mrg_sub_dcdata{ii},sw_sub_dcdata{ii}] ));
    if(tmp_sub_max>ysub_max)
        ysub_max=tmp_sub_max;
    end    
end

dw=0.75;
temp_base=(-dw/2):dw/(nodiam-1):(dw/2);
w_base=repmat(temp_base,[1,nopat]);
w_base=w_base(:);
temp_trans=repmat(1:nodiam,[nopat,1]);
w_trans=temp_trans(:);
xloc_diam=w_base+w_trans; % locations

%% plot

lw=3;
fs=26;
xlab_diam=repmat(diamval,[1,nopat]); % labels
xoff=0.03;
atrans = 0.9;

subplot(2,1,1);
hold on;
for ii=1:nopat
    
    min_epi_mrgth=min( mrg_epi_dcdata{ii} );
    max_epi_mrgth=max( mrg_epi_dcdata{ii} );
    min_epi_swth=min( sw_epi_dcdata{ii} );
    max_epi_swth=max( sw_epi_dcdata{ii} );
    
    % model hresholds
    for jj=1:nodiam
        kk=nodiam*(ii-1)+jj;
        plot([xloc_diam(kk)-xoff,xloc_diam(kk)-xoff],...
            [min_epi_mrgth(jj),max_epi_mrgth(jj)],'k-','linewidth',lw);
        plot([xloc_diam(kk)+xoff,xloc_diam(kk)+xoff],...
            [min_epi_swth(jj),max_epi_swth(jj)],'k--','linewidth',lw);        
    end
    
    % clinical thresholds
    xa=xloc_diam(nodiam*(ii-1)+1)-xoff;
    xb=xloc_diam(nodiam*ii)+xoff;
    ya=isens_expepi(ii);
    yb=idisc_expepi(ii);    
    h_epi_ii=patch([xa,xb,xb,xa],...
                   [ya,ya,yb,yb],'g');
    alpha(h_epi_ii,atrans);

end
hold off;
set(gca,'yscale','log','xtick',[],'ytick',10.^(-1:3),'fontsize',fs);
xlim([xloc_diam(1)-4*xoff,xloc_diam(end)+2*xoff]);
ylim([yepi_min,yepi_max]);

subplot(2,1,2);
hold on;
for ii=1:nopat
    
    min_sub_mrgth=min( mrg_sub_dcdata{ii} );
    max_sub_mrgth=max( mrg_sub_dcdata{ii} );
    min_sub_swth=min( sw_sub_dcdata{ii} );
    max_sub_swth=max( sw_sub_dcdata{ii} );
    
    % model thresholds
    for jj=1:nodiam
        kk=nodiam*(ii-1)+jj;
        plot([xloc_diam(kk)-xoff,xloc_diam(kk)-xoff],...
            [min_sub_mrgth(jj),max_sub_mrgth(jj)],'k-','linewidth',lw);
        plot([xloc_diam(kk)+xoff,xloc_diam(kk)+xoff],...
            [min_sub_swth(jj),max_sub_swth(jj)],'k--','linewidth',lw);        
    end
    
    % clinical thresholds
    xa=xloc_diam(nodiam*(ii-1)+1)-xoff;
    xb=xloc_diam(nodiam*ii)+xoff;
    ya=isens_expsub(ii);
    yb=idisc_expsub(ii);
    h_sub_ii=patch([xa,xb,xb,xa],...
                   [ya,ya,yb,yb],'g');
    alpha(h_sub_ii,atrans);    
    
end
hold off;
set(gca,'yscale','log','ytick',10.^(-1:2));
set(gca,'xtick',xloc_diam,'xticklabel',xlab_diam,'fontsize',fs);
xlim([xloc_diam(1)-4*xoff,xloc_diam(end)+2*xoff]);
ylim([ysub_min,ysub_max]);